!
!     PROGRAMME DIFSOL
!
!     Differences entre les solutions de la libration de la Lune :
!     Sol.LLIB245 - Sol.LLIB04
!     Sol.LLIB403 - Sol.LLIB04
!     Sol.LLIB405 - Sol.LLIB04
!
! --- Declarations -----------------------------------------------------
!
      implicit real*8 (a-h,o-z)
!
      character*50  fname
      dimension     v0(3),v1(3),v2(3),v3(3),dv(3,3)
!
! --- Parametres du calcul ---------------------------------------------
!
      dj0=2360400.5d0
      nbj=109000
      pas=20.d0
!
! --- Fichiers ---------------------------------------------------------
!
      nul1=11
      fname='LLIB04.DAT'
      open (unit=nul1,file=fname,status='old')
      nul2=12
      fname='LLIB245.DAT'
      open (unit=nul2,file=fname,status='old')
      nul3=13
      fname='LLIB403.DAT'
      open (unit=nul3,file=fname,status='old')
      nul4=14
      fname='LLIB405.DAT'
      open (unit=nul4,file=fname,status='old')
      nulout=20
      fname='DIFSOL.TXT'
      open (nulout,file=fname)
!
! --- Calcul des differences -------------------------------------------
!
      npas=nbj/pas
      dj=dj0
!
      do n=1,npas
         if (mod(n,10).eq.0) 
     &   write (*,'(2x,i5.5," / ",i5.5,2x,f9.1)') n,npas,dj
         call LLIB04  (dj,nul1,v0)
         call LLIB245 (dj,nul2,v1)
         call LLIB403 (dj,nul3,v2)
         call LLIB405 (dj,nul4,v3)
         do i=1,3
            dv(i,1)=v1(i)-v0(i)
            dv(i,2)=v2(i)-v0(i)
            dv(i,3)=v3(i)-v0(i)
         enddo
         dv(3,2)=dv(3,2)+2.2d0
         dv(3,3)=dv(3,3)+3.9d0
         an=(dj-2451545.d0)/365.25d0+2000
         write (unit=nulout,'(f10.3,9f12.5)') an,dv
         dj=dj+pas
      enddo
!
      write (*,*) ' FINI '
!
      stop
      end
!
!
!
      subroutine LLIB04 (tjd,nul,var)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Object -----------------------------------------------------------
!
!     Computation of the lunar libration angles (solution LLIB04).
!
! --- Input ------------------------------------------------------------
!
!     tjd      Julian Date TDB (real*8).
!     nul      Logical unit number of the file LLIB04.DAT (integer).
!              This file has to be opened before the first call.
!
! --- Output -----------------------------------------------------------
!
!     var(3)   Libration angles : p1, p2, tau, in arcsecond (real*8).
!              Selenodetic system: axes of principal moments of inertia.              
!
! --- Declarations -----------------------------------------------------
!
      implicit real*8 (a-h,o-z)
!
      real*8,   parameter            :: cpi=3.141592653589793d0
      real*8,   parameter            :: rad=648000.d0/cpi
      real*8,   parameter            :: deg=cpi/180.d0
!
      logical                        :: xinit=.true.
!
      real*8,   dimension(3)         :: var
!
      real*8,   dimension(5)         :: t
      real*8,   dimension(2,500,3)   :: coef
      integer,  dimension(17,500,3)  :: iarg
      integer,  dimension(17)        :: iw
      integer,  dimension(3)         :: nterm
!
      real*8,   dimension(3,5)       :: w
      real*8,   dimension(5)         :: eart, peri
      real*8,   dimension(8,2)       :: p
      real*8,   dimension(4,5)       :: del
      real*8,   dimension(2)         :: zeta
      real*8,   dimension(3)         :: fli
!
      save
!
! --- Function : conversion between (deg,min,sec) and radians ----------
!
      DMS(ideg,imin,sec)=(ideg+imin/60.d0+sec/3600.d0)*deg
!
! --- Time TDT (unit: Julian century) ----------------------------------
!
      t(1) = 1.d0
      t(2) = (tjd-2451545.0d0)/36525.d0
      t(3) = t(2)*t(2)
      t(4) = t(3)*t(2)
      t(5) = t(4)*t(2)
!
! --- First call : initilization of lunar and earth-moon parameters ----
!
      if (xinit) then
!
         w(1,1)  = DMS(218,18,59.8782d0)  ! Mean mean longitude
         w(1,2)  = 1732559343.3328d0/rad  ! of the moon
         w(1,3)  =         -6.8700d0/rad
         w(1,4)  =   0.6604d-2/rad
         w(1,5)  =  -0.3169d-4/rad
!
         w(2,1)  = DMS( 83,21,11.6518d0)  ! Mean longitude 
         w(2,2)  =   14643420.3304d0/rad  ! of the lunar perigee
         w(2,3)  =        -38.2639d0/rad
         w(2,4)  =  -0.45047d-1/rad
         w(2,5)  =   0.21301d-3/rad
!
         w(3,1)  = DMS(125, 2,40.3265d0)  ! Mean longitude 
         w(3,2)  =   -6967919.8851d0/rad  ! of the lunar ascending node
         w(3,3)  =          6.3593d0/rad
         w(3,4)  =   0.7625d-2/rad
         w(3,5)  =  -0.3586d-4/rad
!
         eart(1) = DMS(100,27,59.1880d0)  ! Mean heliocentric
         eart(2) =  129597742.3016d0/rad  ! mean longitude 
         eart(3) =         -0.0202d0/rad  ! of the earth-moon barcenter
         eart(4) =   0.9d-5/rad
         eart(5) =   0.15d-6/rad
!
         peri(1) = DMS(102,56,14.4136d0)  ! Mean longitude
         peri(2) =       1161.2283d0/rad  ! of the perihelion
         peri(3) =          0.5327d0/rad  ! of the earth-moon barcenter
         peri(4) =  -0.138d-3/rad
         peri(5) =   0.d0
!
         fli(1)  =       44820417.d0/rad  ! Free libration
         fli(2)  =        1736493.d0/rad  ! arguments frequencies
         fli(3)  =    -5364715.227d0/rad
!
         p(1,1)  =   DMS(252,15, 3.25986d0)  ! Planetary longitudes J2000
         p(2,1)  =   DMS(181,58,47.28305d0)
         p(3,1)  =   eart(1)
         p(4,1)  =   DMS(355,25,59.78866d0)
         p(5,1)  =   DMS( 34,21, 5.34212d0)
         p(6,1)  =   DMS( 50, 4,38.89694d0)
         p(7,1)  =   DMS(314, 3,18.01841d0)
         p(8,1)  =   DMS(304,20,55.19575d0)
!
         p(1,2)  =   538101628.68898d0/rad  ! Planetary mean motions
         p(2,2)  =   210664136.43355d0/rad
         p(3,2)  =    eart(2)
         p(4,2)  =    68905077.59284d0/rad
         p(5,2)  =    10925660.42861d0/rad
         p(6,2)  =     4399609.65932d0/rad
         p(7,2)  =     1542481.19393d0/rad
         p(8,2)  =      786550.32074d0/rad
!
         dprec   =  0.d0
         preces  = (5029.0966d0+dprec)/rad  ! Constant of the precession
!
         do i=1,5
            del(4,i) = w(1,i)  - eart(i)    ! Arguments of Delaunay
            del(3,i) = w(1,i)  - w(3,i)
            del(1,i) = w(1,i)  - w(2,i)
            del(2,i) = eart(i) - peri(i)
         enddo
!
         del(4,1) = del(4,1)+cpi
         zeta(1)  = w(1,1)
         zeta(2)  = w(1,2)+preces
!
      endif
!
! --- Fist call : reading the lunar libration series LLIB04 ------------
!
      if (xinit) then
!
         read (nul,'(1x)')
         do iv=1,3
            read (nul,'(10x,i10)') nt
            nterm(iv)=nt
            do n=1,nt
               read (nul,'(5x,2d20.13,16i3,i5)') xs,xc,iw
               coef(1,n,iv)=xs; coef(2,n,iv)=xc
               do i=1,17
                  iarg(i,n,iv)=iw(i)
               enddo
            enddo
         enddo
!
         close (nul)
         xinit=.false.
!
      endif
!
! --- Computation of the libration angles p1, p2 and tau ---------------
!
      do iv=1,3
         nt=nterm(iv)
         var(iv)=0.d0
         do n=1,nt
            xs=coef(1,n,iv) ; xc=coef(2,n,iv)
            do i=1,17
               iw(i)=iarg(i,n,iv)
            enddo
            y=0.d0
            itest=iw(17)
            if (itest == 1 .or. itest == 3) then
               do i=1,4
                  do k=1,5
                     y=y+iw(i)*del(i,k)*t(k)
                  enddo
               enddo
            else
               do i=1,4
                  y=y+iw(i)*(del(i,1)+del(i,2)*t(2))
               enddo
            endif
            do i=5,12
               y=y+iw(i)*(p(i-4,1)+p(i-4,2)*t(2))
            enddo
            y=y+iw(13)*(zeta(1)+zeta(2)*t(2))
            do i=14,16
               y=y+iw(i)*fli(i-13)*t(2)
            enddo
            var(iv)=var(iv)+xs*sin(y)+xc*cos(y)
         enddo
      enddo
!
! --- Additive Poisson terms for tau (10-5-98) -------------------------
!
      y      = 18*p(2,1)-16*p(3,1)-del(1,1)+114.56550d0*deg
     &       + (18*p(2,2)-16*p(3,2)-del(1,2))*t(2)
      var(3) = var(3)+0.25d0*sin(y)*t(2)
      y      = del(2,1)+t(2)*del(2,2)
      var(3) = var(3)-0.23d0*sin(y)*t(2)
!
      return
      end subroutine LLIB04
!
!
!
      subroutine LLIB245 (tjd,nul,var)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Object -----------------------------------------------------------
!
!     Computation of the lunar libration angles (solution LLIB245).
!
! --- Input ------------------------------------------------------------
!
!     tjd      Julian Date TDB (real*8).
!     nul      Logical unit number of the file LLIB245.DAT (integer).
!              This file has to be opened before the first call.
!
! --- Output -----------------------------------------------------------
!
!     var(3)   Libration angles : p1, p2, tau, in arcsecond (real*8).
!              Selenodetic system: axes of principal moments of inertia.              
!
! --- Declarations -----------------------------------------------------
!
      implicit real*8 (a-h,o-z)
!
      real*8,   parameter            :: cpi=3.141592653589793d0
      real*8,   parameter            :: rad=648000.d0/cpi
      real*8,   parameter            :: deg=cpi/180.d0
!
      logical                        :: xinit=.true.
!
      real*8,   dimension(3)         :: var
!
      real*8,   dimension(5)         :: t
      real*8,   dimension(2,500,3)   :: coef
      integer,  dimension(17,500,3)  :: iarg
      integer,  dimension(17)        :: iw
      integer,  dimension(3)         :: nterm
!
      real*8,   dimension(3,5)       :: w
      real*8,   dimension(5)         :: eart, peri
      real*8,   dimension(8,2)       :: p
      real*8,   dimension(4,5)       :: del
      real*8,   dimension(2)         :: zeta
      real*8,   dimension(3)         :: fli
!
      save
!
! --- Function : conversion between (deg,min,sec) and radians ----------
!
      DMS(ideg,imin,sec)=(ideg+imin/60.d0+sec/3600.d0)*deg
!
! --- Time TDT (unit: Julian century) ----------------------------------
!
      t(1) = 1.d0
      t(2) = (tjd-2451545.0d0)/36525.d0
      t(3) = t(2)*t(2)
      t(4) = t(3)*t(2)
      t(5) = t(4)*t(2)
!
! --- First call : initilization of lunar and earth-moon parameters ----
!
      if (xinit) then
!
         w(1,1)  =  DMS(218,18,59.83482d0)  ! Mean mean longitude
         w(1,2)  =  1732559343.35614d0/rad  ! of the moon
         w(1,3)  =           -6.7996d0/rad
         w(1,4)  =           0.6604d-2/rad
         w(1,5)  =          -0.3169d-4/rad
!
         w(2,1)  =  DMS( 83,21,11.60739d0)  ! Mean longitude 
         w(2,2)  =     14643420.3398d0/rad  ! of the lunar perigee
         w(2,3)  =          -38.2639d0/rad
         w(2,4)  =         -0.45047d-1/rad
         w(2,5)  =          0.21301d-3/rad
!
         w(3,1)  =  DMS(125, 2,40.29866d0)  ! Mean longitude 
         w(3,2)  =     -6967919.5396d0/rad  ! of the lunar ascending node
         w(3,3)  =            6.3593d0/rad
         w(3,4)  =           0.7625d-2/rad
         w(3,5)  =          -0.3586d-4/rad
!
         eart(1) =  DMS(100,27,59.14346d0)  ! Mean heliocentric
         eart(2) =    129597742.3078d0/rad  ! mean longitude 
         eart(3) =           -0.0202d0/rad  ! of the earth-moon barcenter
         eart(4) =              0.9d-5/rad
         eart(5) =             0.15d-6/rad
!
         peri(1) =  DMS(102,56,14.37127d0)  ! Mean longitude
         peri(2) =         1161.2283d0/rad  ! of the perihelion
         peri(3) =            0.5327d0/rad  ! of the earth-moon barcenter
         peri(4) =           -0.138d-3/rad
         peri(5) =                    0.d0
!
         fli(1)  =       44820553.89d0/rad  ! Free libration
         fli(2)  =       1736494.752d0/rad  ! arguments frequencies
         fli(3)  =      -5364715.227d0/rad
!
         p(1,1)  =  DMS(252,15, 3.25986d0)  ! Planetary longitudes J2000
         p(2,1)  =  DMS(181,58,47.28305d0)
         p(3,1)  =  eart(1)
         p(4,1)  =  DMS(355,25,59.78866d0)
         p(5,1)  =  DMS( 34,21, 5.34212d0)
         p(6,1)  =  DMS( 50, 4,38.89694d0)
         p(7,1)  =  DMS(314, 3,18.01841d0)
         p(8,1)  =  DMS(304,20,55.19575d0)
!
         p(1,2)  =   538101628.68898d0/rad  ! Planetary mean motions
         p(2,2)  =   210664136.43355d0/rad
         p(3,2)  =   eart(2)
         p(4,2)  =    68905077.59284d0/rad
         p(5,2)  =    10925660.42861d0/rad
         p(6,2)  =     4399609.65932d0/rad
         p(7,2)  =     1542481.19393d0/rad
         p(8,2)  =      786550.32074d0/rad
!
         dprec   =  0.d0
         preces  = (5029.0966d0+dprec)/rad  ! Constant of the precession
!
         do i=1,5
            del(4,i) = w(1,i)  - eart(i)    ! Arguments of Delaunay
            del(3,i) = w(1,i)  - w(3,i)
            del(1,i) = w(1,i)  - w(2,i)
            del(2,i) = eart(i) - peri(i)
         enddo
!
         del(4,1) = del(4,1)+cpi
         zeta(1)  = w(1,1)
         zeta(2)  = w(1,2)+preces
!
      endif
!
! --- Fist call : reading the lunar libration series LLIB245 -----------
!
      if (xinit) then
!
         read (nul,'(1x)')
         do iv=1,3
            read (nul,'(10x,i10)') nt
            nterm(iv)=nt
            do n=1,nt
               read (nul,'(5x,2d20.13,16i3,i5)') xs,xc,iw
               coef(1,n,iv)=xs; coef(2,n,iv)=xc
               do i=1,17
                  iarg(i,n,iv)=iw(i)
               enddo
            enddo
         enddo
!
         close (nul)
         xinit=.false.
!
      endif
!
! --- Computation of the libration angles p1, p2 and tau ---------------
!
      do iv=1,3
         nt=nterm(iv)
         var(iv)=0.d0
         do n=1,nt
            xs=coef(1,n,iv) ; xc=coef(2,n,iv)
            do i=1,17
               iw(i)=iarg(i,n,iv)
            enddo
            y=0.d0
            itest=iw(17)
            if (itest == 1 .or. itest == 3) then
               do i=1,4
                  do k=1,5
                     y=y+iw(i)*del(i,k)*t(k)
                  enddo
               enddo
            else
               do i=1,4
                  y=y+iw(i)*(del(i,1)+del(i,2)*t(2))
               enddo
            endif
            do i=5,12
               y=y+iw(i)*(p(i-4,1)+p(i-4,2)*t(2))
            enddo
            y=y+iw(13)*(zeta(1)+zeta(2)*t(2))
            do i=14,16
               y=y+iw(i)*fli(i-13)*t(2)
            enddo
            var(iv)=var(iv)+xs*sin(y)+xc*cos(y)
         enddo
      enddo
!
! --- Additive Poisson terms for tau (10-5-98) -------------------------
!
      y      = 18*p(2,1)-16*p(3,1)-del(1,1)+114.56550d0*deg
     &       + (18*p(2,2)-16*p(3,2)-del(1,2))*t(2)
      var(3) = var(3)+0.25d0*sin(y)*t(2)
      y      = del(2,1)+t(2)*del(2,2)
      var(3) = var(3)-0.23d0*sin(y)*t(2)
!
      return
      end subroutine LLIB245
!
!
!
      subroutine LLIB403 (tjd,nul,var)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Object -----------------------------------------------------------
!
!     Computation of the lunar libration angles (solution LLIB403).
!
! --- Input ------------------------------------------------------------
!
!     tjd      Julian Date TDB (real*8).
!     nul      Logical unit number of the file LLIB403.DAT (integer).
!              This file has to be opened before the first call.
!
! --- Output -----------------------------------------------------------
!
!     var(3)   Libration angles : p1, p2, tau, in arcsecond (real*8).
!              Selenodetic system: axes of principal moments of inertia.              
!
! --- Declarations -----------------------------------------------------
!
      implicit real*8 (a-h,o-z)
!
      real*8,   parameter            :: cpi=3.141592653589793d0
      real*8,   parameter            :: rad=648000.d0/cpi
      real*8,   parameter            :: deg=cpi/180.d0
!
      logical                        :: xinit=.true.
!
      real*8,   dimension(3)         :: var
!
      real*8,   dimension(5)         :: t
      real*8,   dimension(2,500,3)   :: coef
      integer,  dimension(17,500,3)  :: iarg
      integer,  dimension(17)        :: iw
      integer,  dimension(3)         :: nterm
!
      real*8,   dimension(3,5)       :: w
      real*8,   dimension(5)         :: eart, peri
      real*8,   dimension(8,2)       :: p
      real*8,   dimension(4,5)       :: del
      real*8,   dimension(2)         :: zeta
      real*8,   dimension(3)         :: fli
!
      save
!
! --- Function : conversion between (deg,min,sec) and radians ----------
!
      DMS(ideg,imin,sec)=(ideg+imin/60.d0+sec/3600.d0)*deg
!
! --- Time TDT (unit: Julian century) ----------------------------------
!
      t(1) = 1.d0
      t(2) = (tjd-2451545.0d0)/36525.d0
      t(3) = t(2)*t(2)
      t(4) = t(3)*t(2)
      t(5) = t(4)*t(2)
!
! --- First call : initilization of lunar and earth-moon parameters ----
!
      if (xinit) then
!
         w(1,1)  =  DMS(218,18,59.87484d0)  ! Mean mean longitude
         w(1,2)  =  1732559343.35624d0/rad  ! of the moon
         w(1,3)  =           -6.7772d0/rad
         w(1,4)  =           0.6604d-2/rad
         w(1,5)  =          -0.3169d-4/rad
!
         w(2,1)  =  DMS( 83,21,11.64741d0)  ! Mean longitude 
         w(2,2)  =     14643420.3339d0/rad  ! of the lunar perigee
         w(2,3)  =          -38.2639d0/rad
         w(2,4)  =         -0.45047d-1/rad
         w(2,5)  =          0.21301d-3/rad
!
         w(3,1)  =  DMS(125, 2,40.33737d0)  ! Mean longitude 
         w(3,2)  =     -6967919.5455d0/rad  ! of the lunar ascending node
         w(3,3)  =            6.3593d0/rad
         w(3,4)  =           0.7625d-2/rad
         w(3,5)  =          -0.3586d-4/rad
!
         eart(1) =  DMS(100,27,59.18353d0)  ! Mean heliocentric
         eart(2) =    129597742.3027d0/rad  ! mean longitude 
         eart(3) =           -0.0202d0/rad  ! of the earth-moon barcenter
         eart(4) =              0.9d-5/rad
         eart(5) =             0.15d-6/rad
!
         peri(1) =  DMS(102,56,14.41102d0)  ! Mean longitude
         peri(2) =         1161.2283d0/rad  ! of the perihelion
         peri(3) =            0.5327d0/rad  ! of the earth-moon barcenter
         peri(4) =           -0.138d-3/rad
         peri(5) =                    0.d0
!
         fli(1)  =       44820553.89d0/rad  ! Free libration
         fli(2)  =        1736523.38d0/rad  ! arguments frequencies
         fli(3)  =      -5364708.34d0/rad
!
         p(1,1)  =  DMS(252,15, 3.25986d0)  ! Planetary longitudes J2000
         p(2,1)  =  DMS(181,58,47.28305d0)
         p(3,1)  =  eart(1)
         p(4,1)  =  DMS(355,25,59.78866d0)
         p(5,1)  =  DMS( 34,21, 5.34212d0)
         p(6,1)  =  DMS( 50, 4,38.89694d0)
         p(7,1)  =  DMS(314, 3,18.01841d0)
         p(8,1)  =  DMS(304,20,55.19575d0)
!
         p(1,2)  =   538101628.68898d0/rad  ! Planetary mean motions
         p(2,2)  =   210664136.43355d0/rad
         p(3,2)  =   eart(2)
         p(4,2)  =    68905077.59284d0/rad
         p(5,2)  =    10925660.42861d0/rad
         p(6,2)  =     4399609.65932d0/rad
         p(7,2)  =     1542481.19393d0/rad
         p(8,2)  =      786550.32074d0/rad
!
         dprec   =  0.d0
         preces  = (5029.0966d0+dprec)/rad  ! Constant of the precession
!
         do i=1,5
            del(4,i) = w(1,i)  - eart(i)    ! Arguments of Delaunay
            del(3,i) = w(1,i)  - w(3,i)
            del(1,i) = w(1,i)  - w(2,i)
            del(2,i) = eart(i) - peri(i)
         enddo
!
         del(4,1) = del(4,1)+cpi
         zeta(1)  = w(1,1)
         zeta(2)  = w(1,2)+preces
!
      endif
!
! --- Fist call : reading the lunar libration series LLIB403 -----------
!
      if (xinit) then
!
         read (nul,'(1x)')
         do iv=1,3
            read (nul,'(10x,i10)') nt
            nterm(iv)=nt
            do n=1,nt
               read (nul,'(5x,2d20.13,16i3,i5)') xs,xc,iw
               coef(1,n,iv)=xs; coef(2,n,iv)=xc
               do i=1,17
                  iarg(i,n,iv)=iw(i)
               enddo
            enddo
         enddo
!
         close (nul)
         xinit=.false.
!
      endif
!
! --- Computation of the libration angles p1, p2 and tau ---------------
!
      do iv=1,3
         nt=nterm(iv)
         var(iv)=0.d0
         do n=1,nt
            xs=coef(1,n,iv) ; xc=coef(2,n,iv)
            do i=1,17
               iw(i)=iarg(i,n,iv)
            enddo
            y=0.d0
            itest=iw(17)
            if (itest == 1 .or. itest == 3) then
               do i=1,4
                  do k=1,5
                     y=y+iw(i)*del(i,k)*t(k)
                  enddo
               enddo
            else
               do i=1,4
                  y=y+iw(i)*(del(i,1)+del(i,2)*t(2))
               enddo
            endif
            do i=5,12
               y=y+iw(i)*(p(i-4,1)+p(i-4,2)*t(2))
            enddo
            y=y+iw(13)*(zeta(1)+zeta(2)*t(2))
            do i=14,16
               y=y+iw(i)*fli(i-13)*t(2)
            enddo
            var(iv)=var(iv)+xs*sin(y)+xc*cos(y)
         enddo
      enddo
!
! --- Additive Poisson terms for tau (10-5-98) -------------------------
!
      y      = 18*p(2,1)-16*p(3,1)-del(1,1)+114.56550d0*deg
     &       + (18*p(2,2)-16*p(3,2)-del(1,2))*t(2)
      var(3) = var(3)+0.25d0*sin(y)*t(2)
      y      = del(2,1)+t(2)*del(2,2)
      var(3) = var(3)-0.23d0*sin(y)*t(2)
!
      return
      end subroutine LLIB403
!
!
!
      subroutine LLIB405 (tjd,nul,var)
!-----------------------------------------------------------------------
!
!     Ref : JCGF 0407
!
! --- Object -----------------------------------------------------------
!
!     Computation of the lunar libration angles (solution LLIB405).
!
! --- Input ------------------------------------------------------------
!
!     tjd      Julian Date TDB (real*8).
!     nul      Logical unit number of the file LLIB405.DAT (integer).
!              This file has to be opened before the first call.
!
! --- Output -----------------------------------------------------------
!
!     var(3)   Libration angles : p1, p2, tau, in arcsecond (real*8).
!              Selenodetic system: axes of principal moments of inertia.              
!
! --- Declarations -----------------------------------------------------
!
      implicit real*8 (a-h,o-z)
!
      real*8,   parameter            :: cpi=3.141592653589793d0
      real*8,   parameter            :: rad=648000.d0/cpi
      real*8,   parameter            :: deg=cpi/180.d0
!
      logical                        :: xinit=.true.
!
      real*8,   dimension(3)         :: var
!
      real*8,   dimension(5)         :: t
      real*8,   dimension(2,500,3)   :: coef
      integer,  dimension(17,500,3)  :: iarg
      integer,  dimension(17)        :: iw
      integer,  dimension(3)         :: nterm
!
      real*8,   dimension(3,5)       :: w
      real*8,   dimension(5)         :: eart, peri
      real*8,   dimension(8,2)       :: p
      real*8,   dimension(4,5)       :: del
      real*8,   dimension(2)         :: zeta
      real*8,   dimension(3)         :: fli
!
      save
!
! --- Function : conversion between (deg,min,sec) and radians ----------
!
      DMS(ideg,imin,sec)=(ideg+imin/60.d0+sec/3600.d0)*deg
!
! --- Time TDT (unit: Julian century) ----------------------------------
!
      t(1) = 1.d0
      t(2) = (tjd-2451545.0d0)/36525.d0
      t(3) = t(2)*t(2)
      t(4) = t(3)*t(2)
      t(5) = t(4)*t(2)
!
! --- First call : initilization of lunar and earth-moon parameters ----
!
      if (xinit) then
!
         w(1,1)  =  DMS(218,18,59.87267d0)  ! Mean mean longitude
         w(1,2)  =  1732559343.32953d0/rad  ! of the moon
         w(1,3)  =           -6.8368d0/rad
         w(1,4)  =           0.6604d-2/rad
         w(1,5)  =          -0.3169d-4/rad
!
         w(2,1)  =  DMS( 83,21,11.65181d0)  ! Mean longitude 
         w(2,2)  =     14643420.3443d0/rad  ! of the lunar perigee
         w(2,3)  =          -38.2631d0/rad
         w(2,4)  =         -0.45047d-1/rad
         w(2,5)  =          0.21301d-3/rad
!
         w(3,1)  =  DMS(125, 2,40.33472d0)  ! Mean longitude 
         w(3,2)  =     -6967919.5451d0/rad  ! of the lunar ascending node
         w(3,3)  =            6.3590d0/rad
         w(3,4)  =           0.7625d-2/rad
         w(3,5)  =          -0.3586d-4/rad
!
         eart(1) =  DMS(100,27,59.18557d0)  ! Mean heliocentric
         eart(2) =    129597742.3020d0/rad  ! mean longitude 
         eart(3) =           -0.0202d0/rad  ! of the earth-moon barcenter
         eart(4) =              0.9d-5/rad
         eart(5) =             0.15d-6/rad
!
         peri(1) =  DMS(102,56,14.44680d0)  ! Mean longitude
         peri(2) =         1161.2283d0/rad  ! of the perihelion
         peri(3) =            0.5327d0/rad  ! of the earth-moon barcenter
         peri(4) =           -0.138d-3/rad
         peri(5) =                    0.d0
!
         fli(1)  =       44819747.03d0/rad  ! Free libration
         fli(2)  =        1736520.37d0/rad  ! arguments frequencies
         fli(3)  =      -5364600.06d0/rad
!
         p(1,1)  =  DMS(252,15, 3.25986d0)  ! Planetary longitudes J2000
         p(2,1)  =  DMS(181,58,47.28305d0)
         p(3,1)  =  eart(1)
         p(4,1)  =  DMS(355,25,59.78866d0)
         p(5,1)  =  DMS( 34,21, 5.34212d0)
         p(6,1)  =  DMS( 50, 4,38.89694d0)
         p(7,1)  =  DMS(314, 3,18.01841d0)
         p(8,1)  =  DMS(304,20,55.19575d0)
!
         p(1,2)  =   538101628.68898d0/rad  ! Planetary mean motions
         p(2,2)  =   210664136.43355d0/rad
         p(3,2)  =   eart(2)
         p(4,2)  =    68905077.59284d0/rad
         p(5,2)  =    10925660.42861d0/rad
         p(6,2)  =     4399609.65932d0/rad
         p(7,2)  =     1542481.19393d0/rad
         p(8,2)  =      786550.32074d0/rad
!
         dprec   =  0.d0
         preces  = (5029.0966d0+dprec)/rad  ! Constant of the precession
!
         do i=1,5
            del(4,i) = w(1,i)  - eart(i)    ! Arguments of Delaunay
            del(3,i) = w(1,i)  - w(3,i)
            del(1,i) = w(1,i)  - w(2,i)
            del(2,i) = eart(i) - peri(i)
         enddo
!
         del(4,1) = del(4,1)+cpi
         zeta(1)  = w(1,1)
         zeta(2)  = w(1,2)+preces
!
      endif
!
! --- Fist call : reading the lunar libration series LLIB405 -----------
!
      if (xinit) then
!
         read (nul,'(1x)')
         do iv=1,3
            read (nul,'(10x,i10)') nt
            nterm(iv)=nt
            do n=1,nt
               read (nul,'(5x,2d20.13,16i3,i5)') xs,xc,iw
               coef(1,n,iv)=xs; coef(2,n,iv)=xc
               do i=1,17
                  iarg(i,n,iv)=iw(i)
               enddo
            enddo
         enddo
!
         close (nul)
         xinit=.false.
!
      endif
!
! --- Computation of the libration angles p1, p2 and tau ---------------
!
      do iv=1,3
         nt=nterm(iv)
         var(iv)=0.d0
         do n=1,nt
            xs=coef(1,n,iv) ; xc=coef(2,n,iv)
            do i=1,17
               iw(i)=iarg(i,n,iv)
            enddo
            y=0.d0
            itest=iw(17)
            if (itest == 1 .or. itest == 3) then
               do i=1,4
                  do k=1,5
                     y=y+iw(i)*del(i,k)*t(k)
                  enddo
               enddo
            else
               do i=1,4
                  y=y+iw(i)*(del(i,1)+del(i,2)*t(2))
               enddo
            endif
            do i=5,12
               y=y+iw(i)*(p(i-4,1)+p(i-4,2)*t(2))
            enddo
            y=y+iw(13)*(zeta(1)+zeta(2)*t(2))
            do i=14,16
               y=y+iw(i)*fli(i-13)*t(2)
            enddo
            var(iv)=var(iv)+xs*sin(y)+xc*cos(y)
         enddo
      enddo
!
! --- Additive Poisson terms for tau (10-5-98) -------------------------
!
      y      = 18*p(2,1)-16*p(3,1)-del(1,1)+114.56550d0*deg
     &       + (18*p(2,2)-16*p(3,2)-del(1,2))*t(2)
      var(3) = var(3)+0.25d0*sin(y)*t(2)
      y      = del(2,1)+t(2)*del(2,2)
      var(3) = var(3)-0.23d0*sin(y)*t(2)
!
      return
      end subroutine LLIB405
